%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Wave packet propagation integrating the Schrödinger equation,
%  2D time-splitting spectral method,
%  evanescent wave potential.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Wavepacket

%set(1,'DoubleBuffer','on');		mov=avifile('C:\Dokumente und Einstellungen\Philippe\Eigene Dateien\ex3.avi','Quality',100);
addpath C:\Users\Philippe\Physik\Appendix\Appendix;		Constant;
a  =-20e-6;							b  =20e-6;									%grid limits
M  =128;							dz =(b-a)/M;								%grid points
j  =[0:M-1];						Z  =a+j*dz;									%discretization of position space
[ x, z]=meshgrid(Z,Z);
l  =[-M/2:M/2-1];					K  =2*pi*l/(b-a);							%discretization of momentum space
[kx,kz]=meshgrid(K,K);

V  =hP*2*pi* 500*exp(z/5e-6);
E  =hP*2*pi*1000;																		%particle energy
Kx =sqrt(2*87*uA*E)/hP*cos(80*pi/180);											%wavevector
Kz =sqrt(2*87*uA*E)/hP*sin(80*pi/180);
Psi=exp(-((x+15e-6)/4e-6).^2-((z+0e-6)/4e-6).^2).*exp(i*Kx*x+i*Kz*z);					%initial wavefunction
N  =2e5;							g1D=hP*0e-3;

t  =0;								dt =50e-6;									%time
skp=1;
for kk=1:200
	t     =t+dt;
	Psi_x =exp(-i/hP*(V+N*g1D*abs(Psi).^2)*dt/2) .* Psi;
	Psi_xf= fft2(Psi_x);
	Psi_y =exp(i*hP/2/87/uA*(kx.^2+kz.^2)*dt) .* Psi_xf;
	Psi_yf=ifft2(Psi_y);
	Psi   =exp(-i/hP*(V+N*g1D*abs(Psi).^2)*dt/2) .* Psi_yf;
	if((kk/skp)==floor(kk/skp))
		Spx(kk/skp)=sum(Z.*sum(abs(Psi).^2,1) )/sum(sum(abs(Psi).^2,1));		%wavepacket center of mass
		Spz(kk/skp)=sum(Z.*sum(abs(Psi).^2,2)')/sum(sum(abs(Psi).^2,2));
		Spxx(kk/skp)=sqrt(sum(Z.^2.*sum(abs(Psi).^2,1) )/sum(sum(abs(Psi).^2,1)));	%wavepacket width
		Spzz(kk/skp)=sqrt(sum(Z.^2.*sum(abs(Psi).^2,2)')/sum(sum(abs(Psi).^2,2)));

		figure(1);					clf;
		Wavepacket_Plot(t,Z,V/2/pi/hP*1e-3  + abs(Psi).^2*20 + E/2/pi/hP*1e-3,Spx,Spz,1);
		Wavepacket_Quiv(Z,dz,Psi,Spx(kk/skp),Spz(kk/skp));
		Wavepacket_Disp(dt,Spx,Spz,Spxx,Spzz);

%		F  =getframe(gcf);			mov=addframe(mov,F);
		pause(.01);
	end
end
%mov = close(mov);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting subroutine

function Wavepacket_Plot(t,Z,what,Spx,Spz,where)

subplot(2,2,where);					hold on;
surf(Z*1e6,Z*1e6,what);
plot3(Spx*1e6,Spz*1e6,20+0*Spx,'.y','MarkerSize',2);
hold off;
set(gca,'xLim',[min(Z) max(Z)]*1e6,'yLim',[min(Z) max(Z)]*1e6,'FontSize',16,'LineWidth',1.4);
xlabel('\fontsize{18}{\itx}   (\mum)');
ylabel('\fontsize{18}{\itz}   (\mum)');
zlabel('\fontsize{18}{\itV/h}   (kHz)');
if where==1 title(['\fontsize{12}{\itt} = ',num2str(t*1e3),' ms']); end
if where==2 title('\fontsize{12}evanescent matter wave'); end
caxis([0 22]);						shading flat;
view([0 90]);						axis square;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gradient subroutine

function Wavepacket_Quiv(Z,dz,Psi,Spx,Spz)

Constant;
[GPsix,GPsiz]=gradient(     Psi ,dz);
[CPsix,CPsiz]=gradient(conj(Psi),dz);
Flx=-i*hP/2/87/uA*(conj(Psi).*GPsix-Psi.*CPsix);
Flz=-i*hP/2/87/uA*(conj(Psi).*GPsiz-Psi.*CPsiz);

rng=5e-6;
Xx =Z(  (Z>(Spx-rng)) & (Z<(Spx+rng)));
Zz =Z(  (Z>(Spz-rng)) & (Z<(Spz+rng)));
Flx=Flx((Z>(Spz-rng)) & (Z<(Spz+rng)) , (Z>(Spx-rng)) & (Z<(Spx+rng)));
Flz=Flz((Z>(Spz-rng)) & (Z<(Spz+rng)) , (Z>(Spx-rng)) & (Z<(Spx+rng)));
Xx =Xx( 1:2:length(Flz));
Zz =Zz( 1:2:length(Flz));
Flx=Flx(1:2:length(Flz),1:2:length(Flz));
Flz=Flz(1:2:length(Flz),1:2:length(Flz));

subplot(223);						quiver(Xx*1e6,Zz*1e6,Flx*1e6,Flz*1e6,2);
set(gca,'xLim',[Spx-rng Spx+rng]*1e6,'yLim',[Spz-rng Spz+rng]*1e6,'FontSize',16,'LineWidth',1.4);
xlabel('\fontsize{18}{\itx}   (\mum)');
ylabel('\fontsize{18}{\itz}   (\mum)');
axis square;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dispersion subroutine

function Wavepacket_Disp(dt,Spx,Spz,Spxx,Spzz)

tt =[1:length(Spxx)]/length(Spxx)*dt;
subplot(2,2,4);						plot(tt*1e6,Spx*1e6,'r',tt*1e6,Spz*1e6,'g',tt*1e6,Spxx*1e6,'-.r',tt*1e6,Spzz*1e6,'-.g');
set(gca,'FontSize',16,'LineWidth',1.4);
xlabel('\fontsize{18}{\itt}   (\mus)');
ylabel('\fontsize{18}{\itx, z}, \Delta{\itx}, \Delta{\itz}   (\mum)');
